Differential Abundance Analysis (DAA): Fecal microbiome
Eva Kim
11/20/2024
Based on the alpha and beta diversity results, I am not expecting to find a huge difference between groups in terms of bacterial composition. Please note that many of these code are from the Microbiome Intensive Course (MIC Course, Duke University, 2023) and thus the work of Dr. Pixu Shi.
This is after talking with Dr. Shi about what analyses I want to run.
1 Set-Up
First, update the ANCOMBC package if necessary.
#if (!requireNamespace("BiocManager", quietly=TRUE))
# install.packages("BiocManager")
#BiocManager::install("ANCOMBC")2 Load libraries and directories
# rm(list = ls(all = TRUE)) # Unload all packages if necessary
library(phyloseq)
library(fs)
library(rstatix)
library(ggpubr)
library(microViz)
library(ANCOMBC)
library(pals)
library(RColorBrewer)
library(colorBlindness)
library(microbiome)
library(ggrepel)
library(DT)
library(plotly)
library(DESeq2)
library(parallel)
library(dplyr)
library(tidyr)
library(readr)
library(tibble)
library(stringr)rm(list=ls()) # clear environment
git.dir = file.path("/hpc/group/dmc/Eva/GLP_KO_Microbiome_final") # CHANGE ME
ps.rds <- file.path(git.dir, "ps.GLP1.rds")
ps <- read_rds(ps.rds)
ps## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 3304 taxa and 234 samples ]
## sample_data() Sample Data: [ 234 samples by 42 sample variables ]
## tax_table() Taxonomy Table: [ 3304 taxa by 6 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 3304 tips and 3303 internal nodes ]
## refseq() DNAStringSet: [ 3304 reference sequences ]fig.dir = file.path(git.dir, "figures")
map.file = file.path(git.dir, "Metadata_GLP1.csv")
meta.df = read_csv(map.file, show_col_types = FALSE)3 Prepare phyloseq object
I will now aggregate at the genus level, as it is the smallest taxonomic rank that analyses will be run on. In other words, I do not plan on running differential abundance analysis at ASV level.
Let’s prepare the phyloseq object for differential abundance analysis (DAA). Since we are interested in 1) differentially abundant taxa and 2) correlation analysis, we will limit the taxa to only ones that appear in > 3 samples for each sample type.
To avoid running analyses on extremely rare taxa that may be sequencing artifacts, the threshold for minimum number of counts is set at 10. This is, of course, arbitrary. 10 was chosen as 1% of number of read threshold of 1000 above for pruning samples. Based on how I wrote the code, taxa must be present have > 10 reads or more in > 3 samples.
Furthermore, categorical variables should be re-leveled as factors.
# Subset to true fecal samples only and remove tree
ps %>% subset_samples(Sample_type == "Fecal" & Type == "True Sample") -> ps_fecal
ps_fecal_notree <- phyloseq(otu_table(ps_fecal), tax_table(ps_fecal), sample_data(ps_fecal))
# Categorical variables should be factors
ps_fecal_notree@sam_data$Genotype <- factor(ps_fecal_notree@sam_data$Genotype, levels = c("WT", "KO"))
ps_fecal_notree@sam_data$Sex <- factor(ps_fecal_notree@sam_data$Sex, levels = c("Female", "Male"))
ps_fecal_notree@sam_data$Intranasal_Treatment <- factor(ps_fecal_notree@sam_data$Intranasal_Treatment, levels = c("PBS", "HDM"))
ps_fecal_notree@sam_data$Surgery <- factor(ps_fecal_notree@sam_data$Surgery, levels = c("None", "Sham", "VSG"))
ps_fecal_notree@sam_data$Mouse %>% unique() -> mouse_levels
ps_fecal_notree@sam_data$Mouse <- factor(ps_fecal_notree@sam_data$Mouse, levels = mouse_levels)
ps_fecal_notree@sam_data$Diet <- factor(ps_fecal_notree@sam_data$Diet, levels = c("Normal_Chow", "High_Fat_Diet"))
ps_fecal_notree@sam_data$Group <- factor(ps_fecal_notree@sam_data$Group, levels = c("NA", "Control", "1", "2"))
ps_fecal_notree@sam_data$Week <- factor(ps_fecal_notree@sam_data$Week, levels = c("0", "10", "13"))
ps_fecal_notree## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 3304 taxa and 178 samples ]
## sample_data() Sample Data: [ 178 samples by 42 sample variables ]
## tax_table() Taxonomy Table: [ 3304 taxa by 6 taxonomic ranks ]# Aggregate at genus level
ps_g_notree_fecal <- tax_glom(ps_fecal_notree, "Genus")
taxa_names(ps_g_notree_fecal) <- tax_table(ps_g_notree_fecal)[, 'Genus']
ps_g_notree_fecal## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 314 taxa and 178 samples ]
## sample_data() Sample Data: [ 178 samples by 42 sample variables ]
## tax_table() Taxonomy Table: [ 314 taxa by 6 taxonomic ranks ]4 Getting ready: prune data
We will prune the phyloseq object at each phylogenetic rank at and above genus. Since we plan on running correlation analysis after differential abundance, let’s test for taxa that are present in at least 5 samples (arbitrary number) with greater than 10 reads in each sample (also arbitrary. It’s also 1% of 1000 reads, which was the threshold used to remove samples from analysis).
x` Genus
# Pruning: set threshold number
threshold_num = 4
nreads_prune = 10
# Prune taxa at each phylogenetic rank
ps_g_notree_fecal %>%
filter_taxa(., function(x) {sum(x > nreads_prune) > threshold_num}, prune = TRUE) -> ps_genus.prune
ps_genus.prune## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 70 taxa and 178 samples ]
## sample_data() Sample Data: [ 178 samples by 42 sample variables ]
## tax_table() Taxonomy Table: [ 70 taxa by 6 taxonomic ranks ]ps_g_notree_fecal %>% transform_sample_counts(function(x) x/sum(x)) -> ps_genus.ts
ps_genus.ts## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 314 taxa and 178 samples ]
## sample_data() Sample Data: [ 178 samples by 42 sample variables ]
## tax_table() Taxonomy Table: [ 314 taxa by 6 taxonomic ranks ]subset_taxa(ps_genus.ts, ps_genus.ts@tax_table[, "Genus"] %in% ps_genus.prune@tax_table[, "Genus"]) -> ps_genus.prune.ts
ps_genus.prune.ts## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 70 taxa and 178 samples ]
## sample_data() Sample Data: [ 178 samples by 42 sample variables ]
## tax_table() Taxonomy Table: [ 70 taxa by 6 taxonomic ranks ]rowSums(ps_genus.prune.ts@otu_table) %>% as.data.frame() %>%
dplyr::rename(RelAb = ".") %>% mutate(RelAb = RelAb * 100,
RelAb = format(RelAb, scientific = FALSE),
RelAb = as.numeric(RelAb)) %>%
dplyr::arrange(RelAb) -> genus_prune_RelAb
summary(genus_prune_RelAb$RelAb)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 93.13 99.92 99.97 99.89 100.00 100.004.1 Family
ps_fecal_notree %>% tax_glom("Family") -> ps_Family
taxa_names(ps_Family) <- tax_table(ps_Family)[, 'Family']
ps_Family## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 169 taxa and 178 samples ]
## sample_data() Sample Data: [ 178 samples by 42 sample variables ]
## tax_table() Taxonomy Table: [ 169 taxa by 6 taxonomic ranks ]ps_Family %>% transform_sample_counts(function(x) x/sum(x)) -> ps_Family.ts
ps_Family.ts## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 169 taxa and 178 samples ]
## sample_data() Sample Data: [ 178 samples by 42 sample variables ]
## tax_table() Taxonomy Table: [ 169 taxa by 6 taxonomic ranks ]# Prune taxa at each phylogenetic rank
ps_Family %>%
filter_taxa(., function(x) {sum(x > nreads_prune) > threshold_num}, prune = TRUE) -> ps_Family.prune
ps_Family.prune## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 37 taxa and 178 samples ]
## sample_data() Sample Data: [ 178 samples by 42 sample variables ]
## tax_table() Taxonomy Table: [ 37 taxa by 6 taxonomic ranks ]subset_taxa(ps_Family.ts, ps_Family.ts@tax_table[, "Family"] %in% ps_Family.prune@tax_table[, "Family"]) -> ps_Family.prune.ts
ps_Family.prune.ts## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 37 taxa and 178 samples ]
## sample_data() Sample Data: [ 178 samples by 42 sample variables ]
## tax_table() Taxonomy Table: [ 37 taxa by 6 taxonomic ranks ]rowSums(ps_Family.prune.ts@otu_table) %>% as.data.frame() %>%
dplyr::rename(RelAb = ".") %>% mutate(RelAb = RelAb * 100,
RelAb = format(RelAb, scientific = FALSE),
RelAb = as.numeric(RelAb)) %>%
dplyr::arrange(RelAb) -> Family_prune_RelAb
summary(Family_prune_RelAb)## RelAb
## Min. : 94.21
## 1st Qu.: 99.99
## Median :100.00
## Mean : 99.95
## 3rd Qu.:100.00
## Max. :100.004.2 Order
ps_fecal_notree %>% tax_glom("Order") -> ps_Order
taxa_names(ps_Order) <- tax_table(ps_Order)[, 'Order']
ps_Order## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 109 taxa and 178 samples ]
## sample_data() Sample Data: [ 178 samples by 42 sample variables ]
## tax_table() Taxonomy Table: [ 109 taxa by 6 taxonomic ranks ]ps_Order %>% transform_sample_counts(function(x) x/sum(x)) -> ps_Order.ts
ps_Order.ts## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 109 taxa and 178 samples ]
## sample_data() Sample Data: [ 178 samples by 42 sample variables ]
## tax_table() Taxonomy Table: [ 109 taxa by 6 taxonomic ranks ]# Prune taxa at each phylogenetic rank
ps_Order %>%
filter_taxa(., function(x) {sum(x > nreads_prune) > threshold_num}, prune = TRUE) -> ps_Order.prune
ps_Order.prune## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 26 taxa and 178 samples ]
## sample_data() Sample Data: [ 178 samples by 42 sample variables ]
## tax_table() Taxonomy Table: [ 26 taxa by 6 taxonomic ranks ]subset_taxa(ps_Order.ts, ps_Order.ts@tax_table[, "Order"] %in% ps_Order.prune@tax_table[, "Order"]) -> ps_Order.prune.ts
ps_Order.prune.ts## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 26 taxa and 178 samples ]
## sample_data() Sample Data: [ 178 samples by 42 sample variables ]
## tax_table() Taxonomy Table: [ 26 taxa by 6 taxonomic ranks ]rowSums(ps_Order.prune.ts@otu_table) %>% as.data.frame() %>%
dplyr::rename(RelAb = ".") %>% mutate(RelAb = RelAb * 100,
RelAb = format(RelAb, scientific = FALSE),
RelAb = as.numeric(RelAb)) %>%
dplyr::arrange(RelAb) -> Order_prune_RelAb
summary(Order_prune_RelAb)## RelAb
## Min. : 94.22
## 1st Qu.:100.00
## Median :100.00
## Mean : 99.95
## 3rd Qu.:100.00
## Max. :100.004.3 Class
ps_fecal_notree %>% tax_glom("Class") -> ps_Class
taxa_names(ps_Class) <- tax_table(ps_Class)[, 'Class']
ps_Class## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 42 taxa and 178 samples ]
## sample_data() Sample Data: [ 178 samples by 42 sample variables ]
## tax_table() Taxonomy Table: [ 42 taxa by 6 taxonomic ranks ]ps_Class %>% transform_sample_counts(function(x) x/sum(x)) -> ps_Class.ts
ps_Class.ts## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 42 taxa and 178 samples ]
## sample_data() Sample Data: [ 178 samples by 42 sample variables ]
## tax_table() Taxonomy Table: [ 42 taxa by 6 taxonomic ranks ]# Prune taxa at each phylogenetic rank
ps_Class %>%
filter_taxa(., function(x) {sum(x > nreads_prune) > threshold_num}, prune = TRUE) -> ps_Class.prune
ps_Class.prune## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 12 taxa and 178 samples ]
## sample_data() Sample Data: [ 178 samples by 42 sample variables ]
## tax_table() Taxonomy Table: [ 12 taxa by 6 taxonomic ranks ]subset_taxa(ps_Class.ts, ps_Class.ts@tax_table[, "Class"] %in% ps_Class.prune@tax_table[, "Class"]) -> ps_Class.prune.ts
ps_Class.prune.ts## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 12 taxa and 178 samples ]
## sample_data() Sample Data: [ 178 samples by 42 sample variables ]
## tax_table() Taxonomy Table: [ 12 taxa by 6 taxonomic ranks ]rowSums(ps_Class.prune.ts@otu_table) %>% as.data.frame() %>%
dplyr::rename(RelAb = ".") %>% mutate(RelAb = RelAb * 100,
RelAb = format(RelAb, scientific = FALSE),
RelAb = as.numeric(RelAb)) %>%
dplyr::arrange(RelAb) -> Class_prune_RelAb
summary(Class_prune_RelAb)## RelAb
## Min. : 99.99
## 1st Qu.:100.00
## Median :100.00
## Mean :100.00
## 3rd Qu.:100.00
## Max. :100.004.4 Phylum
ps_fecal_notree %>% tax_glom("Phylum") -> ps_Phylum
taxa_names(ps_Phylum) <- tax_table(ps_Phylum)[, 'Phylum']
ps_Phylum## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 20 taxa and 178 samples ]
## sample_data() Sample Data: [ 178 samples by 42 sample variables ]
## tax_table() Taxonomy Table: [ 20 taxa by 6 taxonomic ranks ]ps_Phylum %>% transform_sample_counts(function(x) x/sum(x)) -> ps_Phylum.ts
ps_Phylum.ts## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 20 taxa and 178 samples ]
## sample_data() Sample Data: [ 178 samples by 42 sample variables ]
## tax_table() Taxonomy Table: [ 20 taxa by 6 taxonomic ranks ]# Prune taxa at each phylogenetic rank
ps_Phylum %>%
filter_taxa(., function(x) {sum(x > nreads_prune) > threshold_num}, prune = TRUE) -> ps_Phylum.prune
ps_Phylum.prune## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 10 taxa and 178 samples ]
## sample_data() Sample Data: [ 178 samples by 42 sample variables ]
## tax_table() Taxonomy Table: [ 10 taxa by 6 taxonomic ranks ]subset_taxa(ps_Phylum.ts, ps_Phylum.ts@tax_table[, "Phylum"] %in% ps_Phylum.prune@tax_table[, "Phylum"]) -> ps_Phylum.prune.ts
ps_Phylum.prune.ts## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 10 taxa and 178 samples ]
## sample_data() Sample Data: [ 178 samples by 42 sample variables ]
## tax_table() Taxonomy Table: [ 10 taxa by 6 taxonomic ranks ]rowSums(ps_Phylum.prune.ts@otu_table) %>% as.data.frame() %>%
dplyr::rename(RelAb = ".") %>% mutate(RelAb = RelAb * 100,
RelAb = format(RelAb, scientific = FALSE),
RelAb = as.numeric(RelAb)) %>%
dplyr::arrange(RelAb) -> Phylum_prune_RelAb
summary(Phylum_prune_RelAb)## RelAb
## Min. :100
## 1st Qu.:100
## Median :100
## Mean :100
## 3rd Qu.:100
## Max. :1005 Visualization
Let’s visualize at Phylum level to see how much information was kept:
plot_bar(ps_genus.prune.ts, x = "Label", fill="Phylum") +
facet_grid(scales="free", space = "free_x") +
geom_bar(aes(color=Phylum, fill=Phylum), stat="identity", position="stack") +
theme_bw() +
coord_cartesian(ylim = c(0,1), expand=0) +
theme(axis.text.x=element_blank()) + # not showing samples names
labs(y="Relative Abundance", title = "Relative Abundance of fecal Samples, Pruned @Genus", x = "Sample")+
scale_color_manual(values=SteppedSequential5Steps) +
scale_fill_manual(values=SteppedSequential5Steps) 
plot_bar(ps_Family.prune.ts, x = "Label", fill="Phylum") +
facet_grid(scales="free", space = "free_x") +
geom_bar(aes(color=Phylum, fill=Phylum), stat="identity", position="stack") +
theme_bw() +
coord_cartesian(ylim = c(0,1), expand=0) +
theme(axis.text.x=element_blank()) + # not showing samples names
labs(y="Relative Abundance", title = "Relative Abundance of fecal Samples, Pruned @Family", x = "Sample")+
scale_color_manual(values=SteppedSequential5Steps) +
scale_fill_manual(values=SteppedSequential5Steps) 
plot_bar(ps_Order.prune.ts, x = "Label", fill="Phylum") +
facet_grid(scales="free", space = "free_x") +
geom_bar(aes(color=Phylum, fill=Phylum), stat="identity", position="stack") +
theme_bw() +
coord_cartesian(ylim = c(0,1), expand=0) +
theme(axis.text.x=element_blank()) + # not showing samples names
labs(y="Relative Abundance", title = "Relative Abundance of fecal Samples, Pruned @Order", x = "Sample")+
scale_color_manual(values=SteppedSequential5Steps) +
scale_fill_manual(values=SteppedSequential5Steps) 
plot_bar(ps_Class.prune.ts, x = "Label", fill="Class") +
facet_grid(scales="free", space = "free_x") +
geom_bar(aes(color=Class, fill=Class), stat="identity", position="stack") +
theme_bw() +
coord_cartesian(ylim = c(0,1), expand=0) +
theme(axis.text.x=element_blank()) + # not showing samples names
labs(y="Relative Abundance", title = "Relative Abundance of fecal Samples, Pruned @Class", x = "Sample")+
scale_color_manual(values=SteppedSequential5Steps) +
scale_fill_manual(values=SteppedSequential5Steps) 
plot_bar(ps_Phylum.prune.ts, x = "Label", fill="Phylum") +
facet_grid(scales="free", space = "free_x") +
geom_bar(aes(color=Phylum, fill=Phylum), stat="identity", position="stack") +
theme_bw() +
coord_cartesian(ylim = c(0,1), expand=0) +
theme(axis.text.x=element_blank()) + # not showing samples names
labs(y="Relative Abundance", title = "Relative Abundance of fecal Samples, Pruned @Phylum", x = "Sample")+
scale_color_manual(values=SteppedSequential5Steps) +
scale_fill_manual(values=SteppedSequential5Steps) 
6 Subset metadata
meta.df %>% filter(Type == "True Sample" & Sample_type == "Fecal") %>%
filter(Label %in% phyloseq::sample_data(ps_g_notree_fecal)$Label) -> meta.fecal
meta.fecal7 Q1: Fecal microbiome, HFD, wk 10 & 13, Surgery * Week
Let’s see how many mice we have:
meta.fecal %>% filter(Diet == "High_Fat_Diet" & Week != 0) -> meta.f.hfd
meta.f.hfd %>%
dplyr::count(Week, Surgery, Genotype) %>%
arrange(n)meta.f.hfd %>%
dplyr::count(Week, Surgery) %>%
arrange(n)I will run Week * Surgery.
nrow(meta.f.hfd)## [1] 787.1 Make TSE objects
ps_genus.prune %>% subset_samples(Diet == "High_Fat_Diet" & Week !="0") %>%
mia::makeTreeSummarizedExperimentFromPhyloseq() -> tse.fecal.Q1.Genus
tse.fecal.Q1.Genus## class: TreeSummarizedExperiment
## dim: 70 78
## metadata(0):
## assays(1): counts
## rownames(70): Peptococcus Anaeroplasma ... Bifidobacterium
## Corynebacterium
## rowData names(6): Kingdom Phylum ... Family Genus
## colnames(78): SRR29054904 SRR29054909 ... SRR29055192 SRR29055193
## colData names(42): Label Qubit_ng_ul ... Concentration..ng.uL. is.neg
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## rowLinks: NULL
## rowTree: NULL
## colLinks: NULL
## colTree: NULLps_Family.prune %>% subset_samples(Diet == "High_Fat_Diet" & Week !="0") %>%
mia::makeTreeSummarizedExperimentFromPhyloseq() -> tse.fecal.Q1.Family
tse.fecal.Q1.Family## class: TreeSummarizedExperiment
## dim: 37 78
## metadata(0):
## assays(1): counts
## rownames(37): Acholeplasmataceae Peptococcaceae ... Bifidobacteriaceae
## Corynebacteriaceae
## rowData names(6): Kingdom Phylum ... Family Genus
## colnames(78): SRR29054904 SRR29054909 ... SRR29055192 SRR29055193
## colData names(42): Label Qubit_ng_ul ... Concentration..ng.uL. is.neg
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## rowLinks: NULL
## rowTree: NULL
## colLinks: NULL
## colTree: NULLps_Order.prune %>% subset_samples(Diet == "High_Fat_Diet" & Week !="0") %>%
mia::makeTreeSummarizedExperimentFromPhyloseq() -> tse.fecal.Q1.Order
tse.fecal.Q1.Order## class: TreeSummarizedExperiment
## dim: 26 78
## metadata(0):
## assays(1): counts
## rownames(26): Acholeplasmatales Peptococcales ... Bifidobacteriales
## Corynebacteriales
## rowData names(6): Kingdom Phylum ... Family Genus
## colnames(78): SRR29054904 SRR29054909 ... SRR29055192 SRR29055193
## colData names(42): Label Qubit_ng_ul ... Concentration..ng.uL. is.neg
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## rowLinks: NULL
## rowTree: NULL
## colLinks: NULL
## colTree: NULLps_Class.prune %>% subset_samples(Diet == "High_Fat_Diet" & Week !="0") %>%
mia::makeTreeSummarizedExperimentFromPhyloseq() -> tse.fecal.Q1.Class
tse.fecal.Q1.Class## class: TreeSummarizedExperiment
## dim: 12 78
## metadata(0):
## assays(1): counts
## rownames(12): Campylobacteria Desulfovibrionia ... Clostridia
## Actinobacteria
## rowData names(6): Kingdom Phylum ... Family Genus
## colnames(78): SRR29054904 SRR29054909 ... SRR29055192 SRR29055193
## colData names(42): Label Qubit_ng_ul ... Concentration..ng.uL. is.neg
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## rowLinks: NULL
## rowTree: NULL
## colLinks: NULL
## colTree: NULLps_Phylum.prune %>% subset_samples(Diet == "High_Fat_Diet" & Week !="0") %>%
mia::makeTreeSummarizedExperimentFromPhyloseq() -> tse.fecal.Q1.Phylum
tse.fecal.Q1.Phylum## class: TreeSummarizedExperiment
## dim: 10 78
## metadata(0):
## assays(1): counts
## rownames(10): Campylobacterota Desulfobacterota ... Firmicutes
## Actinobacteriota
## rowData names(6): Kingdom Phylum ... Family Genus
## colnames(78): SRR29054904 SRR29054909 ... SRR29055192 SRR29055193
## colData names(42): Label Qubit_ng_ul ... Concentration..ng.uL. is.neg
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## rowLinks: NULL
## rowTree: NULL
## colLinks: NULL
## colTree: NULL7.2 ANCOM2
7.2.1 Genus
set.seed(123) # Set seed for reproducibility
n_cl_requested = 24
# Genus level
fecal.Q1.Genus <- ANCOMBC::ancom(data=tse.fecal.Q1.Genus, assay_name="counts", tax_level="Genus",
p_adj_method="holm", prv_cut=0, lib_cut=0,
main_var="Surgery", adj_formula="Week",
alpha = 0.05, struc_zero=TRUE,
n_cl = n_cl_requested) ## 'ancom' has been fully evolved to 'ancombc2'.
## Explore the enhanced capabilities of our refined method!res.fecal.Q1.Genus = fecal.Q1.Genus$res
res.fecal.Q1.Genus %>% filter(detected_0.7 == TRUE) -> tab.fecal.Q1.Genus.sig
tab.fecal.Q1.Genus.sig 7.3 Family
set.seed(123) # Set seed for reproducibility
# Family level
fecal.Q1.Family <- ANCOMBC::ancom(data=tse.fecal.Q1.Family, assay_name="counts", tax_level="Family",
p_adj_method="holm", prv_cut=0, lib_cut=0,
main_var="Surgery", adj_formula="Week",
alpha = 0.05, struc_zero=TRUE,
n_cl = n_cl_requested) ## 'ancom' has been fully evolved to 'ancombc2'.
## Explore the enhanced capabilities of our refined method!res.fecal.Q1.Family = fecal.Q1.Family$res
res.fecal.Q1.Family %>% filter(detected_0.7 == TRUE) -> tab.fecal.Q1.Family.sig
tab.fecal.Q1.Family.sig 7.3.1 Order
set.seed(123) # Set seed for reproducibility
# Order level
fecal.Q1.Order <- ANCOMBC::ancom(data=tse.fecal.Q1.Order, assay_name="counts", tax_level="Order",
p_adj_method="holm", prv_cut=0, lib_cut=0,
main_var="Surgery", adj_formula="Week",
alpha = 0.05, struc_zero=TRUE,
n_cl = n_cl_requested) ## 'ancom' has been fully evolved to 'ancombc2'.
## Explore the enhanced capabilities of our refined method!res.fecal.Q1.Order = fecal.Q1.Order$res
res.fecal.Q1.Order %>% filter(detected_0.7 == TRUE) -> tab.fecal.Q1.Order.sig
tab.fecal.Q1.Order.sig 7.3.2 Class
set.seed(123) # Set seed for reproducibility
# Class level
fecal.Q1.Class <- ANCOMBC::ancom(data=tse.fecal.Q1.Class, assay_name="counts", tax_level="Class",
p_adj_method="holm", prv_cut=0, lib_cut=0,
main_var="Surgery", adj_formula="Week",
alpha = 0.05, struc_zero=TRUE,
n_cl = n_cl_requested) ## 'ancom' has been fully evolved to 'ancombc2'.
## Explore the enhanced capabilities of our refined method!res.fecal.Q1.Class = fecal.Q1.Class$res
res.fecal.Q1.Class %>% filter(detected_0.7 == TRUE) -> tab.fecal.Q1.Class.sig
tab.fecal.Q1.Class.sig 7.3.3 Phylum
set.seed(123) # Set seed for reproducibility
# Phylum level
fecal.Q1.Phylum <- ANCOMBC::ancom(data=tse.fecal.Q1.Phylum, assay_name="counts", tax_level="Phylum",
p_adj_method="holm", prv_cut=0, lib_cut=0,
main_var="Surgery", adj_formula="Week",
alpha = 0.05, struc_zero=TRUE,
n_cl = n_cl_requested) ## 'ancom' has been fully evolved to 'ancombc2'.
## Explore the enhanced capabilities of our refined method!res.fecal.Q1.Phylum = fecal.Q1.Phylum$res
res.fecal.Q1.Phylum %>% filter(detected_0.7 == TRUE) -> tab.fecal.Q1.Phylum.sig
tab.fecal.Q1.Phylum.sig 7.3.4 Aggregate Results
tab.fecal.Q1.Genus.sig$Tax_rank <- "Genus"
tab.fecal.Q1.Family.sig$Tax_rank <- "Family"
tab.fecal.Q1.Order.sig$Tax_rank <- "Order"
tab.fecal.Q1.Class.sig$Tax_rank <- "Class"
tab.fecal.Q1.Phylum.sig$Tax_rank <- "Phylum"
rbind(tab.fecal.Q1.Genus.sig,
tab.fecal.Q1.Family.sig,
tab.fecal.Q1.Order.sig,
tab.fecal.Q1.Class.sig,
tab.fecal.Q1.Phylum.sig) -> tab.fecal.Q1
tab.fecal.Q1$Test <- "ANCOM II, fecal microbiome: all HFD, wk 10 & 13, Surgery * Week. Structural zeroes from Surgery"
tab.fecal.Q1 7.4 DESeq2
7.4.1 Genus
ds.fecal.Q1_Genus <- DESeq2::DESeqDataSet(tse.fecal.Q1.Genus, ~ Surgery * Week)## converting counts to integer modedds.fecal.Q1_Genus <- DESeq2::DESeq(ds.fecal.Q1_Genus)## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testing## -- replacing outliers and refitting for 8 genes
## -- DESeq argument 'minReplicatesForReplace' = 7
## -- original counts are preserved in counts(dds)## estimating dispersions## fitting model and testingDESeq2::summary(dds.fecal.Q1_Genus)## $samples
## # A tibble: 1 × 6
## total_counts min_counts max_counts median_counts mean_counts stdev_counts
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1015262 2530 52223 12554. 13016. 6735.
##
## $features
## # A tibble: 1 × 3
## total singletons per_sample_avg
## <int> <int> <dbl>
## 1 70 0 39.2resultsNames(dds.fecal.Q1_Genus)## [1] "Intercept" "Surgery_Sham_vs_None" "Surgery_VSG_vs_None"
## [4] "Week_13_vs_10" "SurgerySham.Week13" "SurgeryVSG.Week13"results(dds.fecal.Q1_Genus, tidy = TRUE, contrast = c("Week", "13", "10")) %>% # tidy converts to df
filter(padj < 0.05) -> fecal.Q1.Genus.Res
fecal.Q1.Genus.Res$Comparison <- c("Week: week 13 over week 10")
results(dds.fecal.Q1_Genus, contrast = c("Surgery", "Sham", "None"), tidy = TRUE) %>%
filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Surgery: Sham over None")
rbind(tmp, fecal.Q1.Genus.Res) -> fecal.Q1.Genus.Res
results(dds.fecal.Q1_Genus, contrast = c("Surgery", "VSG", "None"), tidy = TRUE) %>%
filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Surgery: VSG over None")
rbind(tmp, fecal.Q1.Genus.Res) -> fecal.Q1.Genus.Res
results(dds.fecal.Q1_Genus, name = c("SurgerySham.Week13"), tidy = TRUE) %>%
filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Surgery Sham, wk 13")
rbind(tmp, fecal.Q1.Genus.Res) -> fecal.Q1.Genus.Res
results(dds.fecal.Q1_Genus, name = c("SurgeryVSG.Week13"), tidy = TRUE) %>%
filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Surgery VSG, wk 13")
rbind(tmp, fecal.Q1.Genus.Res) -> fecal.Q1.Genus.Res
fecal.Q1.Genus.Res$Compare_taxon <- "Genus"
fecal.Q1.Genus.Res7.4.2 Family
ds.fecal.Q1_Family <- DESeq2::DESeqDataSet(tse.fecal.Q1.Family, ~ Surgery * Week)## converting counts to integer modedds.fecal.Q1_Family <- DESeq2::DESeq(ds.fecal.Q1_Family)## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testing## -- replacing outliers and refitting for 5 genes
## -- DESeq argument 'minReplicatesForReplace' = 7
## -- original counts are preserved in counts(dds)## estimating dispersions## fitting model and testingDESeq2::summary(dds.fecal.Q1_Family)## $samples
## # A tibble: 1 × 6
## total_counts min_counts max_counts median_counts mean_counts stdev_counts
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1280900 3584 68847 15486. 16422. 8868.
##
## $features
## # A tibble: 1 × 3
## total singletons per_sample_avg
## <int> <int> <dbl>
## 1 37 0 22.8resultsNames(dds.fecal.Q1_Family)## [1] "Intercept" "Surgery_Sham_vs_None" "Surgery_VSG_vs_None"
## [4] "Week_13_vs_10" "SurgerySham.Week13" "SurgeryVSG.Week13"results(dds.fecal.Q1_Family, tidy = TRUE, contrast = c("Week", "13", "10")) %>% # tidy converts to df
filter(padj < 0.05) -> fecal.Q1.Family.Res
fecal.Q1.Family.Res$Comparison <- c("Week: week 13 over week 10")
results(dds.fecal.Q1_Family, contrast = c("Surgery", "Sham", "None"), tidy = TRUE) %>%
filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Surgery: Sham over None")
rbind(tmp, fecal.Q1.Family.Res) -> fecal.Q1.Family.Res
results(dds.fecal.Q1_Family, contrast = c("Surgery", "VSG", "None"), tidy = TRUE) %>%
filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Surgery: VSG over None")
rbind(tmp, fecal.Q1.Family.Res) -> fecal.Q1.Family.Res
results(dds.fecal.Q1_Family, name = c("SurgerySham.Week13"), tidy = TRUE) %>%
filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Surgery Sham, wk 13")
rbind(tmp, fecal.Q1.Family.Res) -> fecal.Q1.Family.Res
results(dds.fecal.Q1_Family, name = c("SurgeryVSG.Week13"), tidy = TRUE) %>%
filter(padj < 0.05) fecal.Q1.Family.Res$Compare_taxon <- "Family"
fecal.Q1.Family.Res7.4.3 Order
ds.fecal.Q1_Order <- DESeq2::DESeqDataSet(tse.fecal.Q1.Order, ~ Surgery * Week)## converting counts to integer modedds.fecal.Q1_Order <- DESeq2::DESeq(ds.fecal.Q1_Order)## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testing## -- replacing outliers and refitting for 4 genes
## -- DESeq argument 'minReplicatesForReplace' = 7
## -- original counts are preserved in counts(dds)## estimating dispersions## fitting model and testingDESeq2::summary(dds.fecal.Q1_Order)## $samples
## # A tibble: 1 × 6
## total_counts min_counts max_counts median_counts mean_counts stdev_counts
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1286962 3587 68977 15488. 16500. 8913.
##
## $features
## # A tibble: 1 × 3
## total singletons per_sample_avg
## <int> <int> <dbl>
## 1 26 0 15.2resultsNames(dds.fecal.Q1_Order)## [1] "Intercept" "Surgery_Sham_vs_None" "Surgery_VSG_vs_None"
## [4] "Week_13_vs_10" "SurgerySham.Week13" "SurgeryVSG.Week13"results(dds.fecal.Q1_Order, tidy = TRUE, contrast = c("Week", "13", "10")) %>% # tidy converts to df
filter(padj < 0.05) -> fecal.Q1.Order.Res
fecal.Q1.Order.Res$Comparison <- c("Week: week 13 over week 10")
results(dds.fecal.Q1_Order, contrast = c("Surgery", "Sham", "None"), tidy = TRUE) %>%
filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Surgery: Sham over None")
rbind(tmp, fecal.Q1.Order.Res) -> fecal.Q1.Order.Res
results(dds.fecal.Q1_Order, contrast = c("Surgery", "VSG", "None"), tidy = TRUE) %>%
filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Surgery: VSG over None")
rbind(tmp, fecal.Q1.Order.Res) -> fecal.Q1.Order.Res
results(dds.fecal.Q1_Order, name = c("SurgerySham.Week13"), tidy = TRUE) %>%
filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Surgery Sham, wk 13")
rbind(tmp, fecal.Q1.Order.Res) -> fecal.Q1.Order.Res
results(dds.fecal.Q1_Order, name = c("SurgeryVSG.Week13"), tidy = TRUE) %>%
filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Surgery VSG, wk 13")
rbind(tmp, fecal.Q1.Order.Res) -> fecal.Q1.Order.Res
fecal.Q1.Order.Res$Compare_taxon <- "Order"
fecal.Q1.Order.Res7.4.4 Class
ds.fecal.Q1_Class <- DESeq2::DESeqDataSet(tse.fecal.Q1.Class, ~ Surgery * Week)## converting counts to integer modedds.fecal.Q1_Class <- DESeq2::DESeq(ds.fecal.Q1_Class)## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testing## -- replacing outliers and refitting for 1 genes
## -- DESeq argument 'minReplicatesForReplace' = 7
## -- original counts are preserved in counts(dds)## estimating dispersions## fitting model and testingDESeq2::summary(dds.fecal.Q1_Class)## $samples
## # A tibble: 1 × 6
## total_counts min_counts max_counts median_counts mean_counts stdev_counts
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1288152 3587 69016 15488. 16515. 8912.
##
## $features
## # A tibble: 1 × 3
## total singletons per_sample_avg
## <int> <int> <dbl>
## 1 12 0 9.33resultsNames(dds.fecal.Q1_Class)## [1] "Intercept" "Surgery_Sham_vs_None" "Surgery_VSG_vs_None"
## [4] "Week_13_vs_10" "SurgerySham.Week13" "SurgeryVSG.Week13"results(dds.fecal.Q1_Class, tidy = TRUE, contrast = c("Week", "13", "10")) %>% # tidy converts to df
filter(padj < 0.05) -> fecal.Q1.Class.Res
fecal.Q1.Class.Res$Comparison <- c("Week: week 13 over week 10")
results(dds.fecal.Q1_Class, contrast = c("Surgery", "Sham", "None"), tidy = TRUE) %>%
filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Surgery: Sham over None")
rbind(tmp, fecal.Q1.Class.Res) -> fecal.Q1.Class.Res
results(dds.fecal.Q1_Class, contrast = c("Surgery", "VSG", "None"), tidy = TRUE) %>%
filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Surgery: VSG over None")
rbind(tmp, fecal.Q1.Class.Res) -> fecal.Q1.Class.Res
results(dds.fecal.Q1_Class, name = c("SurgerySham.Week13"), tidy = TRUE) %>%
filter(padj < 0.05) results(dds.fecal.Q1_Class, name = c("SurgeryVSG.Week13"), tidy = TRUE) %>%
filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Surgery VSG, wk 13")
rbind(tmp, fecal.Q1.Class.Res) -> fecal.Q1.Class.Res
fecal.Q1.Class.Res$Compare_taxon <- "Class"
fecal.Q1.Class.Res7.4.5 Phylum
ds.fecal.Q1_Phylum <- DESeq2::DESeqDataSet(tse.fecal.Q1.Phylum, ~ Surgery * Week)## converting counts to integer modedds.fecal.Q1_Phylum <- DESeq2::DESeq(ds.fecal.Q1_Phylum)## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testing## -- replacing outliers and refitting for 1 genes
## -- DESeq argument 'minReplicatesForReplace' = 7
## -- original counts are preserved in counts(dds)## estimating dispersions## fitting model and testingDESeq2::summary(dds.fecal.Q1_Phylum)## $samples
## # A tibble: 1 × 6
## total_counts min_counts max_counts median_counts mean_counts stdev_counts
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1288193 3587 69016 15490. 16515. 8912.
##
## $features
## # A tibble: 1 × 3
## total singletons per_sample_avg
## <int> <int> <dbl>
## 1 10 0 8resultsNames(dds.fecal.Q1_Phylum)## [1] "Intercept" "Surgery_Sham_vs_None" "Surgery_VSG_vs_None"
## [4] "Week_13_vs_10" "SurgerySham.Week13" "SurgeryVSG.Week13"results(dds.fecal.Q1_Phylum, tidy = TRUE, contrast = c("Week", "13", "10")) %>% # tidy converts to df
filter(padj < 0.05) -> fecal.Q1.Phylum.Res
fecal.Q1.Phylum.Res$Comparison <- c("Week: week 13 over week 10")
results(dds.fecal.Q1_Phylum, contrast = c("Surgery", "Sham", "None"), tidy = TRUE) %>%
filter(padj < 0.05) results(dds.fecal.Q1_Phylum, contrast = c("Surgery", "VSG", "None"), tidy = TRUE) %>%
filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Surgery: VSG over None")
rbind(tmp, fecal.Q1.Phylum.Res) -> fecal.Q1.Phylum.Res
results(dds.fecal.Q1_Phylum, name = c("SurgerySham.Week13"), tidy = TRUE) %>%
filter(padj < 0.05) results(dds.fecal.Q1_Phylum, name = c("SurgeryVSG.Week13"), tidy = TRUE) %>%
filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Surgery VSG, wk 13")
rbind(tmp, fecal.Q1.Phylum.Res) -> fecal.Q1.Phylum.Res
fecal.Q1.Phylum.Res$Compare_taxon <- "Phylum"
fecal.Q1.Phylum.Res7.4.6 Aggregate results
rbind(fecal.Q1.Genus.Res,
fecal.Q1.Family.Res,
fecal.Q1.Class.Res,
fecal.Q1.Order.Res,
fecal.Q1.Phylum.Res) -> fecal.Q1.deseq
fecal.Q1.deseq$Test <- "DESeq2, Fecal microbiome: all HFD, wk 10 + 13, Surgery * Week"
fecal.Q1.deseq7.5 Intersection
tab.fecal.Q1 %>% filter(taxon %in% fecal.Q1.deseq$row) -> ancom.intersect.hfd
ancom.intersect.hfdfecal.Q1.deseq %>% filter(row %in% tab.fecal.Q1$taxon) 7.6 Graph
psmelt(ps_genus.prune.ts) -> Genus.prune.melt
psmelt(ps_Family.prune.ts) -> Family.prune.melt
psmelt(ps_Order.prune.ts) -> Order.prune.melt
psmelt(ps_Class.prune.ts) -> Class.prune.melt
psmelt(ps_Phylum.prune.ts) -> Phylum.prune.meltbind_rows(
Genus.prune.melt,
Family.prune.melt,
Order.prune.melt,
Class.prune.melt,
Phylum.prune.melt
) -> fecal.RA.Full
fecal.RA.Full$Surgery <- factor(fecal.RA.Full$Surgery, levels = c("None", "Sham", "VSG"))
fecal.RA.Full %>%
dplyr::select(OTU:Label, Sample_type:Mouse, Genotype:Surgery, Week) -> fecal.RA.Full
head(fecal.RA.Full)tail(fecal.RA.Full)fecal.RA.Full %>% filter(Diet == "High_Fat_Diet" & Week != "0") -> fecal.graph.hfd
head(fecal.graph.hfd)tail(fecal.graph.hfd)surgery.labs <- c("No Surgery", "Sham Surgery", "VSG")
names(surgery.labs) <- c("None", "Sham", "VSG")
fecal.graph.hfd %>%
filter(OTU %in% ancom.intersect.hfd$taxon) %>%
ggplot(aes(x = Surgery, y = Abundance, color = Surgery)) +
geom_boxplot() + geom_point() +
labs(y="Relative Abundance", title = str_glue(ancom.intersect.hfd$taxon[1], " Relative Abundance, fecal microbiome"),
x="Surgery") +
theme_bw(base_size = 14) +
facet_wrap(~Week)
fecal.graph.hfd %>%
filter(OTU %in% ancom.intersect.hfd$taxon) %>%
ggplot(aes(x = Week, y = Abundance, color = Week)) +
geom_boxplot() + geom_point() +
labs(y="Relative Abundance", title = str_glue(ancom.intersect.hfd$taxon[1], " Relative Abundance, fecal microbiome"),
x="Surgery") +
theme_bw(base_size = 14) +
facet_wrap(~Surgery, scales = "free")
fecal.graph.hfd %>%
filter(OTU %in% ancom.intersect.hfd$taxon) %>%
group_by(Week) %>%
kruskal_test(Abundance ~ Surgery) %>%
adjust_pvalue(method = "holm") %>%
add_significance()Post-hoc: pairwise comparison with Wilcoxon test.
fecal.graph.hfd %>%
filter(OTU %in% ancom.intersect.hfd$taxon) %>%
group_by(Week) %>%
wilcox_test(Abundance ~ Surgery, p.adjust.method = "holm") %>%
add_significance() %>%
add_y_position() -> fecal.hfd.res
fecal.hfd.resweek.labs <- c("Week 0", "Week 10", "Week 13")
names(week.labs) <- c("0", "10", "13")
fecal.graph.hfd %>%
filter(OTU %in% ancom.intersect.hfd$taxon) %>%
ggplot(aes(x = Surgery, y = Abundance, color = Surgery)) +
geom_boxplot() + geom_point() +
labs(y="Lactobacillaceae Relative Abundance",
x="Surgery") +
theme_bw(base_size = 14) +
facet_wrap(~Week, labeller = labeller(Week = week.labs)) + theme(legend.position = "top") -> fecal.hfd.res.graph
fecal.hfd.res.graph
fecal.hfd.res$y.position <- c(0.45, 0.5, 0.57, 0.55, 0.64, 0.7)
fecal.hfd.res.graph +
stat_pvalue_manual(fecal.hfd.res, label = "p.adj.signif") -> fecal.hfd.res.graph
fecal.hfd.res.graph
ragg::agg_jpeg(file.path(fig.dir, "fecal_surgery_lactobacillaceae.jpeg"), width = 7, height = 5, units = "in", res = 600)
fecal.hfd.res.graph
dev.off()## png
## 28 Q2: no surgery group, all weeks, Diet * Week + Genotype
Let’s see how many mice we have:
meta.fecal %>% filter(Surgery == "None") -> meta.f.nosurgery
meta.f.nosurgery %>%
dplyr::count(Week, Diet, Genotype) %>%
arrange(n)meta.f.nosurgery %>%
dplyr::count(Week, Diet) %>%
arrange(n)I will run Diet * Week + Genotype.
nrow(meta.f.nosurgery)## [1] 1208.1 Make TSE objects
ps_genus.prune %>% subset_samples(Surgery == "None") %>%
mia::makeTreeSummarizedExperimentFromPhyloseq() -> tse.fecal.Q2.Genus
tse.fecal.Q2.Genus## class: TreeSummarizedExperiment
## dim: 70 120
## metadata(0):
## assays(1): counts
## rownames(70): Peptococcus Anaeroplasma ... Bifidobacterium
## Corynebacterium
## rowData names(6): Kingdom Phylum ... Family Genus
## colnames(120): SRR29054901 SRR29054902 ... SRR29055188 SRR29055191
## colData names(42): Label Qubit_ng_ul ... Concentration..ng.uL. is.neg
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## rowLinks: NULL
## rowTree: NULL
## colLinks: NULL
## colTree: NULLps_Family.prune %>% subset_samples(Surgery == "None") %>%
mia::makeTreeSummarizedExperimentFromPhyloseq() -> tse.fecal.Q2.Family
tse.fecal.Q2.Family## class: TreeSummarizedExperiment
## dim: 37 120
## metadata(0):
## assays(1): counts
## rownames(37): Acholeplasmataceae Peptococcaceae ... Bifidobacteriaceae
## Corynebacteriaceae
## rowData names(6): Kingdom Phylum ... Family Genus
## colnames(120): SRR29054901 SRR29054902 ... SRR29055188 SRR29055191
## colData names(42): Label Qubit_ng_ul ... Concentration..ng.uL. is.neg
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## rowLinks: NULL
## rowTree: NULL
## colLinks: NULL
## colTree: NULLps_Order.prune %>% subset_samples(Surgery == "None") %>%
mia::makeTreeSummarizedExperimentFromPhyloseq() -> tse.fecal.Q2.Order
tse.fecal.Q2.Order## class: TreeSummarizedExperiment
## dim: 26 120
## metadata(0):
## assays(1): counts
## rownames(26): Acholeplasmatales Peptococcales ... Bifidobacteriales
## Corynebacteriales
## rowData names(6): Kingdom Phylum ... Family Genus
## colnames(120): SRR29054901 SRR29054902 ... SRR29055188 SRR29055191
## colData names(42): Label Qubit_ng_ul ... Concentration..ng.uL. is.neg
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## rowLinks: NULL
## rowTree: NULL
## colLinks: NULL
## colTree: NULLps_Class.prune %>% subset_samples(Surgery == "None") %>%
mia::makeTreeSummarizedExperimentFromPhyloseq() -> tse.fecal.Q2.Class
tse.fecal.Q2.Class## class: TreeSummarizedExperiment
## dim: 12 120
## metadata(0):
## assays(1): counts
## rownames(12): Campylobacteria Desulfovibrionia ... Clostridia
## Actinobacteria
## rowData names(6): Kingdom Phylum ... Family Genus
## colnames(120): SRR29054901 SRR29054902 ... SRR29055188 SRR29055191
## colData names(42): Label Qubit_ng_ul ... Concentration..ng.uL. is.neg
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## rowLinks: NULL
## rowTree: NULL
## colLinks: NULL
## colTree: NULLps_Phylum.prune %>% subset_samples(Surgery == "None") %>%
mia::makeTreeSummarizedExperimentFromPhyloseq() -> tse.fecal.Q2.Phylum
tse.fecal.Q2.Phylum## class: TreeSummarizedExperiment
## dim: 10 120
## metadata(0):
## assays(1): counts
## rownames(10): Campylobacterota Desulfobacterota ... Firmicutes
## Actinobacteriota
## rowData names(6): Kingdom Phylum ... Family Genus
## colnames(120): SRR29054901 SRR29054902 ... SRR29055188 SRR29055191
## colData names(42): Label Qubit_ng_ul ... Concentration..ng.uL. is.neg
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## rowLinks: NULL
## rowTree: NULL
## colLinks: NULL
## colTree: NULL8.2 ANCOM II
8.2.1 Genus
set.seed(123) # Set seed for reproducibility
# Genus level
fecal.Q2.Genus <- ANCOMBC::ancom(data=tse.fecal.Q2.Genus, assay_name="counts", tax_level="Genus",
p_adj_method="holm", prv_cut=0, lib_cut=0,
main_var="Diet", adj_formula="Week + Genotype",
alpha = 0.05, struc_zero=TRUE,
n_cl = n_cl_requested) ## 'ancom' has been fully evolved to 'ancombc2'.
## Explore the enhanced capabilities of our refined method!res.fecal.Q2.Genus = fecal.Q2.Genus$res
res.fecal.Q2.Genus %>% filter(detected_0.7 == TRUE) -> tab.fecal.Q2.Genus.sig
tab.fecal.Q2.Genus.sig 8.2.2 Family
set.seed(123) # Set seed for reproducibility
# Family level
fecal.Q2.Family <- ANCOMBC::ancom(data=tse.fecal.Q2.Family, assay_name="counts", tax_level="Family",
p_adj_method="holm", prv_cut=0, lib_cut=0,
main_var="Diet", adj_formula="Week + Genotype",
alpha = 0.05, struc_zero=TRUE,
n_cl = n_cl_requested) ## 'ancom' has been fully evolved to 'ancombc2'.
## Explore the enhanced capabilities of our refined method!res.fecal.Q2.Family = fecal.Q2.Family$res
res.fecal.Q2.Family %>% filter(detected_0.7 == TRUE) -> tab.fecal.Q2.Family.sig
tab.fecal.Q2.Family.sig 8.2.3 Order
set.seed(123) # Set seed for reproducibility
# Order level
fecal.Q2.Order <- ANCOMBC::ancom(data=tse.fecal.Q2.Order, assay_name="counts", tax_level="Order",
p_adj_method="holm", prv_cut=0, lib_cut=0,
main_var="Diet", adj_formula="Week + Genotype",
alpha = 0.05, struc_zero=TRUE,
n_cl = n_cl_requested) ## 'ancom' has been fully evolved to 'ancombc2'.
## Explore the enhanced capabilities of our refined method!res.fecal.Q2.Order = fecal.Q2.Order$res
res.fecal.Q2.Order %>% filter(detected_0.7 == TRUE) -> tab.fecal.Q2.Order.sig
tab.fecal.Q2.Order.sig 8.2.4 Class
set.seed(123) # Set seed for reproducibility
# Class level
fecal.Q2.Class <- ANCOMBC::ancom(data=tse.fecal.Q2.Class, assay_name="counts", tax_level="Class",
p_adj_method="holm", prv_cut=0, lib_cut=0,
main_var="Diet", adj_formula="Week + Genotype",
alpha = 0.05, struc_zero=TRUE,
n_cl = n_cl_requested) ## 'ancom' has been fully evolved to 'ancombc2'.
## Explore the enhanced capabilities of our refined method!res.fecal.Q2.Class = fecal.Q2.Class$res
res.fecal.Q2.Class %>% filter(detected_0.7 == TRUE) -> tab.fecal.Q2.Class.sig
tab.fecal.Q2.Class.sig 8.2.5 Phylum
set.seed(123) # Set seed for reproducibility
# Phylum level
fecal.Q2.Phylum <- ANCOMBC::ancom(data=tse.fecal.Q2.Phylum, assay_name="counts", tax_level="Phylum",
p_adj_method="holm", prv_cut=0, lib_cut=0,
main_var="Diet", adj_formula="Week + Genotype",
alpha = 0.05, struc_zero=TRUE,
n_cl = n_cl_requested) ## 'ancom' has been fully evolved to 'ancombc2'.
## Explore the enhanced capabilities of our refined method!res.fecal.Q2.Phylum = fecal.Q2.Phylum$res
res.fecal.Q2.Phylum %>% filter(detected_0.7 == TRUE) -> tab.fecal.Q2.Phylum.sig
tab.fecal.Q2.Phylum.sig 8.2.6 Aggregate Results
tab.fecal.Q2.Genus.sig$Tax_rank <- "Genus"
tab.fecal.Q2.Family.sig$Tax_rank <- "Family"
tab.fecal.Q2.Order.sig$Tax_rank <- "Order"
tab.fecal.Q2.Class.sig$Tax_rank <- "Class"
tab.fecal.Q2.Phylum.sig$Tax_rank <- "Phylum"
rbind(tab.fecal.Q2.Genus.sig,
tab.fecal.Q2.Family.sig,
tab.fecal.Q2.Order.sig,
tab.fecal.Q2.Class.sig,
tab.fecal.Q2.Phylum.sig) -> tab.fecal.Q2
tab.fecal.Q2$Test <- "ANCOM II, fecal microbiome: no surgery, main effect Diet (structural zeroes) and week + genotype"
tab.fecal.Q2 8.3 DESeq2
8.3.1 Genus
ds.fecal.Q2_Genus <- DESeq2::DESeqDataSet(tse.fecal.Q2.Genus, ~ Diet * Week + Genotype)## converting counts to integer modedds.fecal.Q2_Genus <- DESeq2::DESeq(ds.fecal.Q2_Genus)## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testingDESeq2::summary(dds.fecal.Q2_Genus)## $samples
## # A tibble: 1 × 6
## total_counts min_counts max_counts median_counts mean_counts stdev_counts
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1238163 2007 45853 10035 10318. 5570.
##
## $features
## # A tibble: 1 × 3
## total singletons per_sample_avg
## <int> <int> <dbl>
## 1 70 0 43.2resultsNames(dds.fecal.Q2_Genus)## [1] "Intercept" "Diet_High_Fat_Diet_vs_Normal_Chow"
## [3] "Week_10_vs_0" "Week_13_vs_0"
## [5] "Genotype_KO_vs_WT" "DietHigh_Fat_Diet.Week10"
## [7] "DietHigh_Fat_Diet.Week13"results(dds.fecal.Q2_Genus, tidy = TRUE, name =c("Diet_High_Fat_Diet_vs_Normal_Chow")) %>% # tidy converts to df
filter(padj < 0.05) results(dds.fecal.Q2_Genus, name =c("Genotype_KO_vs_WT"), tidy = TRUE) %>%
filter(padj < 0.05) -> fecal.Q2.Genus.Res
fecal.Q2.Genus.Res$Comparison <- c("Genotype: KO vs. WT")
results(dds.fecal.Q2_Genus, name =c("DietHigh_Fat_Diet.Week10"), tidy = TRUE) %>%
filter(padj < 0.05) -> tmp
tmp$Comparison <- c("HFD, wk 10")
rbind(tmp, fecal.Q2.Genus.Res) -> fecal.Q2.Genus.Res
results(dds.fecal.Q2_Genus, name =c("DietHigh_Fat_Diet.Week13"), tidy = TRUE) %>%
filter(padj < 0.05) -> tmp
tmp$Comparison <- c("HFD, wk 13")
rbind(tmp, fecal.Q2.Genus.Res) -> fecal.Q2.Genus.Res
results(dds.fecal.Q2_Genus, name =c("Week_10_vs_0"), tidy = TRUE) %>%
filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Week: 10 vs. 0")
rbind(tmp, fecal.Q2.Genus.Res) -> fecal.Q2.Genus.Res
results(dds.fecal.Q2_Genus, name =c("Week_13_vs_0"), tidy = TRUE) %>%
filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Week: 13 vs. 0")
rbind(tmp, fecal.Q2.Genus.Res) -> fecal.Q2.Genus.Res
fecal.Q2.Genus.Res$Compare_taxon <- "Genus"
fecal.Q2.Genus.Res8.3.2 Family
ds.fecal.Q2_Family <- DESeq2::DESeqDataSet(tse.fecal.Q2.Family, ~ Diet * Week + Genotype)## converting counts to integer modedds.fecal.Q2_Family <- DESeq2::DESeq(ds.fecal.Q2_Family)## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testingDESeq2::summary(dds.fecal.Q2_Family)## $samples
## # A tibble: 1 × 6
## total_counts min_counts max_counts median_counts mean_counts stdev_counts
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1571691 2593 56081 12346. 13097. 6989.
##
## $features
## # A tibble: 1 × 3
## total singletons per_sample_avg
## <int> <int> <dbl>
## 1 37 0 25.3results(dds.fecal.Q2_Family, tidy = TRUE, name =c("Diet_High_Fat_Diet_vs_Normal_Chow")) %>% # tidy converts to df
filter(padj < 0.05) -> fecal.Q2.Family.Res
fecal.Q2.Family.Res$Comparison <- c("Diet: HFD vs. NC")
results(dds.fecal.Q2_Family, name =c("Genotype_KO_vs_WT"), tidy = TRUE) %>%
filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Genotype: KO vs. WT")
rbind(tmp, fecal.Q2.Family.Res) -> fecal.Q2.Family.Res
results(dds.fecal.Q2_Family, name =c("DietHigh_Fat_Diet.Week10"), tidy = TRUE) %>%
filter(padj < 0.05) -> tmp
tmp$Comparison <- c("HFD, wk 10")
rbind(tmp, fecal.Q2.Family.Res) -> fecal.Q2.Family.Res
results(dds.fecal.Q2_Family, name =c("DietHigh_Fat_Diet.Week13"), tidy = TRUE) %>%
filter(padj < 0.05) -> tmp
tmp$Comparison <- c("HFD, wk 13")
rbind(tmp, fecal.Q2.Family.Res) -> fecal.Q2.Family.Res
results(dds.fecal.Q2_Family, name =c("Week_10_vs_0"), tidy = TRUE) %>%
filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Week: 10 vs. 0")
rbind(tmp, fecal.Q2.Family.Res) -> fecal.Q2.Family.Res
results(dds.fecal.Q2_Family, name =c("Week_13_vs_0"), tidy = TRUE) %>%
filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Week: 13 vs. 0")
rbind(tmp, fecal.Q2.Family.Res) -> fecal.Q2.Family.Res
fecal.Q2.Family.Res$Compare_taxon <- "Family"
fecal.Q2.Family.Res8.3.3 Order
ds.fecal.Q2_Order <- DESeq2::DESeqDataSet(tse.fecal.Q2.Order, ~ Diet * Week + Genotype)## converting counts to integer modedds.fecal.Q2_Order <- DESeq2::DESeq(ds.fecal.Q2_Order)## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testingDESeq2::summary(dds.fecal.Q2_Order)## $samples
## # A tibble: 1 × 6
## total_counts min_counts max_counts median_counts mean_counts stdev_counts
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1593986 2595 57077 12560 13283. 7111.
##
## $features
## # A tibble: 1 × 3
## total singletons per_sample_avg
## <int> <int> <dbl>
## 1 26 0 18.0results(dds.fecal.Q2_Order, tidy = TRUE, name =c("Diet_High_Fat_Diet_vs_Normal_Chow")) %>% # tidy converts to df
filter(padj < 0.05) results(dds.fecal.Q2_Order, name =c("Genotype_KO_vs_WT"), tidy = TRUE) %>%
filter(padj < 0.05) -> fecal.Q2.Order.Res
fecal.Q2.Order.Res$Comparison <- c("Genotype: KO vs. WT")
results(dds.fecal.Q2_Order, name =c("DietHigh_Fat_Diet.Week10"), tidy = TRUE) %>%
filter(padj < 0.05) -> tmp
tmp$Comparison <- c("HFD, wk 10")
rbind(tmp, fecal.Q2.Order.Res) -> fecal.Q2.Order.Res
results(dds.fecal.Q2_Order, name =c("DietHigh_Fat_Diet.Week13"), tidy = TRUE) %>%
filter(padj < 0.05) -> tmp
tmp$Comparison <- c("HFD, wk 13")
rbind(tmp, fecal.Q2.Order.Res) -> fecal.Q2.Order.Res
results(dds.fecal.Q2_Order, name =c("Week_10_vs_0"), tidy = TRUE) %>%
filter(padj < 0.05) results(dds.fecal.Q2_Order, name =c("Week_13_vs_0"), tidy = TRUE) %>%
filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Week: 13 vs. 0")
rbind(tmp, fecal.Q2.Order.Res) -> fecal.Q2.Order.Res
fecal.Q2.Order.Res$Compare_taxon <- "Order"
fecal.Q2.Order.Res8.3.4 Class
ds.fecal.Q2_Class <- DESeq2::DESeqDataSet(tse.fecal.Q2.Class, ~ Diet * Week + Genotype)## converting counts to integer modedds.fecal.Q2_Class <- DESeq2::DESeq(ds.fecal.Q2_Class)## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testingDESeq2::summary(dds.fecal.Q2_Class)## $samples
## # A tibble: 1 × 6
## total_counts min_counts max_counts median_counts mean_counts stdev_counts
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1594382 2595 57079 12562 13287. 7112.
##
## $features
## # A tibble: 1 × 3
## total singletons per_sample_avg
## <int> <int> <dbl>
## 1 12 0 10.6results(dds.fecal.Q2_Class, tidy = TRUE, name =c("Diet_High_Fat_Diet_vs_Normal_Chow")) %>% # tidy converts to df
filter(padj < 0.05) results(dds.fecal.Q2_Class, name =c("Genotype_KO_vs_WT"), tidy = TRUE) %>%
filter(padj < 0.05) -> fecal.Q2.Class.Res
fecal.Q2.Class.Res$Comparison <- c("Genotype: KO vs. WT")
results(dds.fecal.Q2_Class, name =c("DietHigh_Fat_Diet.Week10"), tidy = TRUE) %>%
filter(padj < 0.05) -> tmp
tmp$Comparison <- c("HFD, wk 10")
rbind(tmp, fecal.Q2.Class.Res) -> fecal.Q2.Class.Res
results(dds.fecal.Q2_Class, name =c("DietHigh_Fat_Diet.Week13"), tidy = TRUE) %>%
filter(padj < 0.05) -> tmp
tmp$Comparison <- c("HFD, wk 13")
rbind(tmp, fecal.Q2.Class.Res) -> fecal.Q2.Class.Res
results(dds.fecal.Q2_Class, name =c("Week_10_vs_0"), tidy = TRUE) %>%
filter(padj < 0.05) results(dds.fecal.Q2_Class, name =c("Week_13_vs_0"), tidy = TRUE) %>%
filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Week: 13 vs. 0")
rbind(tmp, fecal.Q2.Class.Res) -> fecal.Q2.Class.Res
fecal.Q2.Class.Res$Compare_taxon <- "Class"
fecal.Q2.Class.Res8.3.5 Phylum
ds.fecal.Q2_Phylum <- DESeq2::DESeqDataSet(tse.fecal.Q2.Phylum, ~ Diet * Week + Genotype)## converting counts to integer modedds.fecal.Q2_Phylum <- DESeq2::DESeq(ds.fecal.Q2_Phylum)## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testingDESeq2::summary(dds.fecal.Q2_Phylum)## $samples
## # A tibble: 1 × 6
## total_counts min_counts max_counts median_counts mean_counts stdev_counts
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1594436 2595 57083 12562 13287. 7112.
##
## $features
## # A tibble: 1 × 3
## total singletons per_sample_avg
## <int> <int> <dbl>
## 1 10 0 8.79results(dds.fecal.Q2_Phylum, tidy = TRUE, name =c("Diet_High_Fat_Diet_vs_Normal_Chow")) %>% # tidy converts to df
filter(padj < 0.05) results(dds.fecal.Q2_Phylum, name =c("Genotype_KO_vs_WT"), tidy = TRUE) %>%
filter(padj < 0.05) -> fecal.Q2.Phylum.Res
fecal.Q2.Phylum.Res$Comparison <- c("Genotype: KO vs. WT")
results(dds.fecal.Q2_Phylum, name =c("DietHigh_Fat_Diet.Week10"), tidy = TRUE) %>%
filter(padj < 0.05) -> tmp
tmp$Comparison <- c("HFD, wk 10")
rbind(tmp, fecal.Q2.Phylum.Res) -> fecal.Q2.Phylum.Res
results(dds.fecal.Q2_Phylum, name =c("DietHigh_Fat_Diet.Week13"), tidy = TRUE) %>%
filter(padj < 0.05) -> tmp
tmp$Comparison <- c("HFD, wk 13")
rbind(tmp, fecal.Q2.Phylum.Res) -> fecal.Q2.Phylum.Res
results(dds.fecal.Q2_Phylum, name =c("Week_10_vs_0"), tidy = TRUE) %>%
filter(padj < 0.05) results(dds.fecal.Q2_Phylum, name =c("Week_13_vs_0"), tidy = TRUE) %>%
filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Week: 13 vs. 0")
rbind(tmp, fecal.Q2.Phylum.Res) -> fecal.Q2.Phylum.Res
fecal.Q2.Phylum.Res$Compare_taxon <- "Phylum"
fecal.Q2.Phylum.Res8.3.6 Aggregate results
rbind(fecal.Q2.Genus.Res,
fecal.Q2.Family.Res,
fecal.Q2.Class.Res,
fecal.Q2.Order.Res,
fecal.Q2.Phylum.Res) -> fecal.Q2.deseq
fecal.Q2.deseq$Test <- "DESeq2, Fecal microbiome: no surgery, ~ Diet * Week + Genotype"
fecal.Q2.deseq8.4 Intersection
tab.fecal.Q2 %>% filter(taxon %in% fecal.Q2.deseq$row) -> ancom.intersect.nosurgery
ancom.intersect.nosurgeryfecal.Q2.deseq %>% filter(row %in% tab.fecal.Q2$taxon) %>%
dplyr::select(row, Comparison, Compare_taxon, Test) %>%
group_by(row) %>%
mutate(n.ind = paste0(1:n())) %>%
pivot_wider(
names_from = "n.ind",
values_from = Comparison
) -> fecal.Q2.deseq.intersect.wide
fecal.Q2.deseq.intersect.wide8.5 Graph
fecal.RA.Full %>% filter(Surgery == "None") -> fecal.graph.nosurgery
head(fecal.graph.nosurgery)tail(fecal.graph.nosurgery)for (i in 1:length(ancom.intersect.nosurgery$taxon)) {
fecal.graph.nosurgery %>%
filter(OTU == ancom.intersect.nosurgery$taxon[i]) %>%
ggplot(aes(x = as.character(Week), y = Abundance, color = Diet)) +
geom_boxplot() + geom_point() +
labs(y="Relative Abundance", title = str_glue(ancom.intersect.nosurgery$taxon[i], " Relative Abundance, fecal"),
x="Weel") +
facet_wrap(~Diet) +
theme_bw(base_size = 14) +
theme(legend.position = "top") -> p
print(p)
}
9 Export csv files
fecal.Q1.deseq %>% filter(row %in% ancom.intersect.hfd$taxon) -> fecal.Q1.deseq.intersect
fecal.Q2.deseq %>% filter(row %in% ancom.intersect.nosurgery$taxon) -> fecal.Q2.deseq.intersectwrite_csv(fecal.RA.Full, file.path(git.dir, "fecal_relative_abundance.csv"), append = FALSE, col_names = TRUE)
write_csv(ancom.intersect.hfd, file.path(git.dir, "fecal_ancom_hfd_intersect.csv"), append = FALSE, col_names = TRUE)
write_csv(ancom.intersect.nosurgery, file.path(git.dir, "fecal_ancom_nosurgery_intersect.csv"), append = FALSE, col_names = TRUE)
write_csv(fecal.Q1.deseq.intersect, file.path(git.dir, "fecal_deseq_hfd_intersect.csv"), append = FALSE, col_names = TRUE)
write_csv(fecal.Q2.deseq.intersect, file.path(git.dir, "fecal_deseq_nosurgery_intersect.csv"), append = FALSE, col_names = TRUE)10 Reproducibility
sessionInfo()## R version 4.3.1 (2023-06-16)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.3 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so; LAPACK version 3.10.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Etc/UTC
## tzcode source: system (glibc)
##
## attached base packages:
## [1] parallel stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] stringr_1.5.0 tibble_3.2.1
## [3] readr_2.1.4 tidyr_1.3.0
## [5] dplyr_1.1.3 DESeq2_1.40.2
## [7] SummarizedExperiment_1.30.2 Biobase_2.60.0
## [9] MatrixGenerics_1.12.3 matrixStats_1.1.0
## [11] GenomicRanges_1.52.1 GenomeInfoDb_1.36.4
## [13] IRanges_2.34.1 S4Vectors_0.38.2
## [15] BiocGenerics_0.46.0 plotly_4.10.3
## [17] DT_0.29 ggrepel_0.9.4
## [19] microbiome_1.22.0 colorBlindness_0.1.9
## [21] RColorBrewer_1.1-3 pals_1.8
## [23] ANCOMBC_2.2.2 microViz_0.11.0
## [25] ggpubr_0.6.0 ggplot2_3.4.4
## [27] rstatix_0.7.2 fs_1.6.3
## [29] phyloseq_1.44.0
##
## loaded via a namespace (and not attached):
## [1] bitops_1.0-7 DirichletMultinomial_1.42.0
## [3] httr_1.4.7 doParallel_1.0.17
## [5] numDeriv_2016.8-1.1 tools_4.3.1
## [7] doRNG_1.8.6 backports_1.4.1
## [9] utf8_1.2.4 R6_2.5.1
## [11] vegan_2.6-4 lazyeval_0.2.2
## [13] mgcv_1.9-0 rhdf5filters_1.12.1
## [15] permute_0.9-7 withr_2.5.1
## [17] gridExtra_2.3 textshaping_0.3.7
## [19] cli_3.6.1 sandwich_3.1-0
## [21] labeling_0.4.3 sass_0.4.7
## [23] mvtnorm_1.2-3 proxy_0.4-27
## [25] systemfonts_1.0.5 yulab.utils_0.1.0
## [27] foreign_0.8-85 dichromat_2.0-0.1
## [29] scater_1.28.0 decontam_1.20.0
## [31] maps_3.4.1 readxl_1.4.3
## [33] rstudioapi_0.15.0 RSQLite_2.3.3
## [35] generics_0.1.3 gridGraphics_0.5-1
## [37] vroom_1.6.4 gtools_3.9.4
## [39] car_3.1-2 Matrix_1.6-1.1
## [41] biomformat_1.28.0 ggbeeswarm_0.7.2
## [43] fansi_1.0.5 DescTools_0.99.52
## [45] DECIPHER_2.28.0 abind_1.4-5
## [47] lifecycle_1.0.3 multcomp_1.4-25
## [49] yaml_2.3.7 carData_3.0-5
## [51] rhdf5_2.44.0 Rtsne_0.16
## [53] grid_4.3.1 blob_1.2.4
## [55] crayon_1.5.2 lattice_0.22-5
## [57] beachmat_2.16.0 cowplot_1.1.1
## [59] mapproj_1.2.11 pillar_1.9.0
## [61] knitr_1.44 boot_1.3-28.1
## [63] gld_2.6.6 codetools_0.2-19
## [65] glue_1.6.2 data.table_1.14.8
## [67] MultiAssayExperiment_1.26.0 vctrs_0.6.4
## [69] treeio_1.24.3 Rdpack_2.6
## [71] cellranger_1.1.0 gtable_0.3.4
## [73] cachem_1.0.8 xfun_0.40
## [75] rbibutils_2.2.16 S4Arrays_1.2.1
## [77] survival_3.5-7 SingleCellExperiment_1.22.0
## [79] iterators_1.0.14 gmp_0.7-2
## [81] TH.data_1.1-2 nlme_3.1-163
## [83] bit64_4.0.5 bslib_0.5.1
## [85] irlba_2.3.5.1 vipor_0.4.5
## [87] rpart_4.1.21 colorspace_2.1-0
## [89] DBI_1.1.3 Hmisc_5.1-1
## [91] nnet_7.3-19 ade4_1.7-22
## [93] Exact_3.2 tidyselect_1.2.0
## [95] bit_4.0.5 compiler_4.3.1
## [97] htmlTable_2.4.2 BiocNeighbors_1.18.0
## [99] expm_0.999-7 DelayedArray_0.26.7
## [101] checkmate_2.3.0 scales_1.2.1
## [103] digest_0.6.33 minqa_1.2.6
## [105] rmarkdown_2.25 XVector_0.40.0
## [107] htmltools_0.5.6.1 pkgconfig_2.0.3
## [109] base64enc_0.1-3 lme4_1.1-34
## [111] sparseMatrixStats_1.12.2 fastmap_1.1.1
## [113] rlang_1.1.1 htmlwidgets_1.6.2
## [115] DelayedMatrixStats_1.22.6 farver_2.1.1
## [117] jquerylib_0.1.4 zoo_1.8-12
## [119] jsonlite_1.8.7 energy_1.7-11
## [121] BiocParallel_1.34.2 BiocSingular_1.16.0
## [123] RCurl_1.98-1.13 magrittr_2.0.3
## [125] Formula_1.2-5 scuttle_1.10.3
## [127] GenomeInfoDbData_1.2.10 Rhdf5lib_1.22.1
## [129] munsell_0.5.0 Rcpp_1.0.11
## [131] ape_5.7-1 viridis_0.6.4
## [133] CVXR_1.0-11 stringi_1.7.12
## [135] rootSolve_1.8.2.4 zlibbioc_1.46.0
## [137] MASS_7.3-60 plyr_1.8.9
## [139] lmom_3.0 Biostrings_2.68.1
## [141] splines_4.3.1 multtest_2.56.0
## [143] hms_1.1.3 locfit_1.5-9.8
## [145] igraph_1.5.1 ggsignif_0.6.4
## [147] rngtools_1.5.2 reshape2_1.4.4
## [149] ScaledMatrix_1.8.1 evaluate_0.22
## [151] tzdb_0.4.0 nloptr_2.0.3
## [153] foreach_1.5.2 purrr_1.0.2
## [155] rsvd_1.0.5 broom_1.0.5
## [157] Rmpfr_0.9-3 e1071_1.7-13
## [159] tidytree_0.4.5 ragg_1.2.6
## [161] viridisLite_0.4.2 class_7.3-22
## [163] gsl_2.1-8 lmerTest_3.1-3
## [165] memoise_2.0.1 beeswarm_0.4.0
## [167] cluster_2.1.4 TreeSummarizedExperiment_2.8.0
## [169] mia_1.8.0